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We study the dynamics of topological defects in the triangular Ising antiferromagnet, a related model 
on the square lattice equivalent to the six-vertex ice model, and the three-state antiferromagnetic 
Potts model on the square lattice. Since each of these models has a height representation in which 
defects are screw dislocations, we expect them to be attracted or repelled with an entropy-driven 
Coulomb force. In each case we show explicitly how this force is felt through local fields. We 
measure the force numerically, both by quenching the system to zero temperature and by measuring 
the motion of vortex pairs. For the three-state Potts model, we calculate both the force and the 
defect mobility, and find reasonable agreement with theory. 
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I. INTRODUCTION 

Topological defects play a major role in mediating 
phase transitions. In this paper, we look at several mod- 
els with discrete degrees of freedom: the triangular Ising 
antiferromagnet, a related model on the square lattice 
equivalent to the six-vertex ice model, and the g = 3 
Potts antiferromagnet on the square lattice. In each 
model we review how to describe defects as vortices or 
screw dislocations in a height representation, and write 
the total charge inside a finite region as a sum around its 
perimeter. 

In each case, we find that the motion of the defects 
is proportional to a local field F. If the states far away 
from a defect are uncorrelated, then this field decreases 
as Bjl-KT where B is the Burgers vector of the defect. We 
therefore argue that defects are governed by a first-order 
Coulomb force. 
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acting like vortices in a viscous regime or dislocations in 
a solid. 

To confirm this, we do two kinds of numerical exper- 
iments. In a quench to zero temperature, we measure 
the density of defects as a function of time, and find log- 
arithmic corrections to p ex 1/i like those seen in the 
XY model. However, a free diffusion experiment sug- 
gests that this occurs whenever a local conservation law 
is used to determine the defects' initial positions and den- 
sity fluctuations, whether or not there are forces between 
them. Therefore, we also attempt to measure the force 
directly by placing a pair of defects r apart, allowing the 
lattice around them to equilibrate, and measuring how 
they would move toward or away from each other if we 
allowed them to. We obtain good agreement with a force 



of the form (Ar) = A/{r + rg), although finite-size ef- 
fects and a partial lack of ergodicity at T = affect our 
results. 

Since the energy of a pair of defects does not depend 
on the distance between them, this force must be entrop- 
ically driven. We review how both the height model and 
heuristic arguments give a term in the entropy propor- 
tional to — Inr, whose gradient then gives a 1/r force. 
Using Park and Widom's calculation of the stiffness of 
the equivalent height model, we predict that A for the 
Potts antiferromagnet should be 3/4. This is in reason- 
able agreement with our results. 



II. THE TRIANGULAR ISING 
ANTIFERROMAGNET 



The Ising antiferromagnet 
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on the triangular lattice can be experimentally realized in 
anhydrous alums with a yavapaiite layered structure |q| 
such as RbFe(S04)2. Its equilibrium behavior was solved 
by Wannier [Q and Houtappel IQ. Its critical point is 
at zero temperature, where it has a nonzero entropy per 
site of 
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Stephenson [p2| showed that its correlations at T 
decay algebraically as r"^/^. 

Since every triangular plaquet of the lattice must con- 
tain at least one frustrated bond, we define a defect as a 
plaquet where all three bonds are frustrated, i.e. where 
all three sites are in the same state. An isolated defect 



+ + 

- + - ^^ - + - 
+ - + - 



- + - + 
+ + +^^+ - + 

- + - + 



+ + 

- + - 

+ + 



+ + 



+ + 



FIG. 1. Defects in the triangular antiferromagnet diffusing 
and annihilating. 
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FIG. 2. Two ways of calculating the defect charge inside a 
region. Equation hi gives (l/6)(— 2 — 2— 1 — 1) — —1, and SiWi 
winds once clockwise around the origin. Both show a single 
negative (A) defect. 



cannot be annealed away with local changes; as figure 
m shows, defects diffuse under single spin-flip dynamics, 
maintaining their orientation on the lattice while revers- 
ing their spins. Thus there are two types of defects, V 
and A. We assign these charges of -t-1 and —1 respec- 
tively, and pairs of opposite type can annihilate when 
they meet. These defects seem to have been first noticed 
by Landau Q . 

We can define the charge within a region as an integral 
going counterclockwise along the bonds b of its perimeter. 
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An example is shown in figure 0. Here t is a unit vec- 
tor that lends an orientation to each bond in the lattice, 
and the expression in the parentheses is 2 for bonds con- 
necting like sites and —1 for unlike ones. This latter 
quantity was used by Blote et al. MM as a height differ- 
ence between neighboring sites to map the triangular an- 
tiferromagnet onto a solid-on-solid model of the (1,1,1) 
corner of a simple cubic crystal, and this was generalized 
by Zeng and Henley pS] to arbitrary spin. Defects then 
correspond to screw dislocations with Burgers vector ±6. 
Another way to define the charge is through a local 
magnetization, whose winding is like a vorticity |15|. We 



FIG. 3. The local magnetization sw around a pair of de- 
fects. This picture was obtained by relaxing a 33 x 33 rhombic 
lattice until only two defects remained, and then holding them 
fixed for 10^ updates per site while we averaged the magneti- 
zation at each site. Finally, sw was averaged over triplets of 
neighboring sites to smooth it. 



can three-color the sites of the lattice with vectors 27r/3 



apart, w 



ifc/3 £pj. fc = 1^ 2, such that moving posi- 



tively along t rotates w counterclockwise by 2tt/3. Then 
Q is the winding number of SiWi around the origin, which 
we can write as a contour integral if we like: 
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This method is also shown in figure 0. In figures and 
we show the time average of sw on a lattice with a 
single pair of defects. The fields around both defects, 
with opposite windings, are clearly visible. 

If we define a field E by rotating t clockwise by 7r/2 and 
multiplying the result by the same quantity (3/2)siSi+i + 
(1/2) as in equation |l|, then E will tend to point towards 
negative defects and away from positive ones. We can 
then rewrite equation n^ as 
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We can now establish a linear relationship between E 
and the expected motion of a defect. As figure || shows, 
once we choose which site of a positive defect to fiip, the 
center of the new defect will be equal to the location of 
that site plus E/{2y/3) (we take the lattice spacing to be 
1) where we measure E on the dashed edge facing the 
site being flipped. Since all three sites are equally likely, 
the average movement is 
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FIG. 4. A scatterplot of the averaged values of sw shown 
in figure pi 
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where the sum is over the three dashed edges surrounding 
the defect. 

Now if we draw a perimeter of radius r around a defect 
of charge q, and if we assume the system is isotropic, 
equation impHes that the average field is 
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pointing away from the defect. Here we assume that 
other defects are sufficiently far away so that they don't 
bias the surrounding lattice. If this is true, then combin- 
ing it with equation gives 
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for the expected motion of another defect of charge q'. 
In other words, we expect the average motion to be gov- 
erned by a Coulomb force, in a viscous regime with con- 
stant mobility. The coefficient VS/Stt is only approxi- 
mate, since we have used a local mean-field approxima- 
tion in which E has no correlations between the dashed 
edges and so (1/3) J^E = E. 

Obviously, this will only hold if the field is sufficiently 
uncorrelated at this distance, and if the presence of an- 
other defect adds linearly to E rather than nonlinearly. 
For instance, if curvature of one defect's field is corre- 
lated with the presence of another, the field lines will 
be bunched between the two defects, and the force will 
fall off more slowly than 1/r. In the sections below, we 
perform various numerical experiments to measure (Ax) , 
and we claim that it does in fact decrease as 1/r. 
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FIG. 5. Defining E and establishing its relationship with 
(Ax). There are two cases where E = 1 away from the defect, 
and averaging them gives (A£) — l/(2\/3) towards the dashed 
edge. When E = 2 towards the defect, the site does not fiip 
since this would increase the energy by creating three defects 
instead of one, and relative to that site the defect moves l/-\/3 
back to its original position. In the case not shown, both spins 
on the dashed edge are +, and the defect annihilates with a 
neighboring one. We can ignore this case for defects far apart. 



III. A RELATED MODEL ON THE SQUARE 

LATTICE: THE NO-STRIPES RULE AND THE 

SIX- VERTEX ICE MODEL 

Consider the following model on the square lattice: 
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S1S4 - S2S3 + S1S2S3S4 



If we like, we can think of this as two Ising ferromagnets, 
one on each sublattice, strongly coupled by the four-point 
interaction like a kind of Ashkin- Teller model [nl. The 
Hamiltonian has the effect of giving plaquets of the form 

, , and their rotations and reversals 

energies of —1, while giving plaquets of the form 

and their rotations an energy of -1-3. We call these last 
plaquets 'stripes,' and they act as defects in the system. 
As figure pi shows, isolated stripes diffuse and annihilate 
in pairs. 

This model is related to the triangular antifcrromag- 
net through the charge-preserving zig-zag transformation 
shown in figure |7|. This preserves energy for the most 
part, with the one discrepancy that a defect-free square 
plaquet can be mapped onto an adjacent pair of positive 
and negative defects in the triangular lattice. Since such 
a pair will quickly annihilate, this should not significantly 
affect the dynamics. 

Nevertheless, this 'no stripes' rule is quite elegant, and 
is worth studying on its own. There are two types of 
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FIG. 6. Defects in the triangular Ising antiferromagnet dif- 
fusing and annihilating. 
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FIG. 7. By flipping half the sites in a zig-zag pattern, we 
can transform the 'no stripes' model into the triangular anti- 
ferromagnet in a charge-preserving way. 



stripes, one which is vertical 
squares and horizontal 
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when on even lattice 



when on odd ones, and 
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the other which is vice versa. If we assign charges -1-1 
and —1 to these respectively, the total charge Q inside 
a region can be written as a sum over bonds b going 
counterclockwise around it, 
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-1 if 6 is vertical 
-1 if 5 is horizontal 



+1 if 6 has an even square on its left 
— 1 if & has an odd square on its left 



(5) 



An example is given in figure y. 

If we define a field E perpendicular to each bond by 
rotating b clockwise or counterclockwise tt/2 according to 
whether the quantity summed in equation g| is positive or 
negative, then E is the same regardless of the orientation 
of b. Moreover, E will tend to point toward negative 
stripes and away from positive ones, and we can rewrite 
equation || as 
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FIG. 8. An example of the sum in equation B. The charge 
enclosed is (1/4) ^ = +1. 
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As shown in figure 0, if our single spin-flip dynamics 
chooses the upper left-hand site of a positive stripe, then 
the stripe will move by (—1/2,-1/2) plus one-half the 
sum of E on the two edges radiating outward from that 
site. When the similar statement is made for the other 
corners, the constants (±1/2, ±1/2) cancel, and the av- 
erage movement of a stripe of charge q is 
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where the sum is over the eight edges radiating out from 
the stripe's plaquet, shown dashed in the figure. By the 
same argument as for the triangular antiferromagnet, we 
argue that there is a first-order force that decreases as 
l/r. Specifically, equation H gives 
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if we assume isotropy, and 
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approximation in which 



in a local mean-field 
{l/^)Y.E = E. 

If we look at what combination of forces can occur 
around a plaquet, figure M shows that E maps the no- 
stripes rule onto the six- vertex ice model, where each site 
of the dual lattice has two incoming and two outgoing 
arrows. This can also be thought of as a loop-covering 
model, in which every vertex is covered by exactly one 
loop S. The entropy per site of this model has been 
calculated exactly by Lieb Ea] and is | In | . 

In this mapping, positive and negative stripes cor- 
respond to positive and negative monopoles, with four 
outgoing or incoming arrows. A single spin- flip corre- 
sponds to reversing the four arrows around a site; this 
preserves the eight-vertex ice rule that each site have an 
even number of incoming arrows. If we only allow these 
moves when no new monopoles will be created, they cause 
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FIG. 9. The average movement of a single stripe is propor- 
tional to the average field on the eight dashed edges. Here we 
show movement resulting from flipping the stripe's upper left 
site. In the lower-right hand case, the stripe does not move 
at all, since flipping this site would increase the energy by 
creating three stripes instead of one. 
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FIG. 11. Defining charge in the no-stripes rule as the wind- 
ing number of sw around the origin. The upper-left site is 
even, and the winding of -1-1 corresponds to a positive defect. 

where each site has a state Si E {1,2, ... ,q} and 6 is the 
Kronecker delta function. The states are often thought 
of as colors, so that the ground state consists of a color- 
ing where no two neighbors have the same color. When 
q — 3, this can be thought of as a discretization of the 
antiferromagnetic XY model, 
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FIG. 10. Using E to map the no-stripes rule to the 
six- vertex ice rule. The plaquets shown are even. Positive and 
negative stripes become positive and negative monopoles. 



monopoles to migrate. Thus we claim that under this 
dynamics, monopoles in the six-vertex ice model also in- 
teract with a Coulomb force. 

Knops points out that the no-stripes rule or the six- 
vertex ice model can also be mapped to a Villain model 
pql . We can then define a vector w on the two sublattices 
as we did for the triangular antiferromagnet. Let u; be 1 
and e'^'/^ on odd and even sites respectively. Then Q is 
the winding number of SiWi around the origin, 

Zm J sw 

as shown in figure |llj, just as for the Ising antiferromag- 
net. Finally, this model can also be mapped to a body- 
centered solid-on-solid model [|6|, in which defects corre- 
spond to screw dislocations of Burgers vector ±4. 



IV. THE Q = 3 SQUARE POTTS 
ANTIFERROMAGNET 

The g-state antiferromagnetic Potts model is a gener- 
alization of the Ising antiferromagnet. 
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where each site has a unit spin pointing in some direction 
6. If the 9i are restricted to three directions 27r/3 apart, 
then equations [g and [^ are equal up to an affine trans- 
formation since cos{9i — 9j) depends only on whether 6i 
and 9j are the same or different. 

We will focus on this three-state model on the square 
lattice. Baxter [Q and Nightingale and Schick p7| ] 
showed that it is critical at T = 0, and that this is its 
only phase transition. Baxter solved the T — case as 
a hard-squares model |y. For recent Monte Carlo stud- 
ies using the Wang-Swendsen-Kotecky cluster algorithm 
Q, see e.g. Ferreira and Sokal Q. 

Kolafa pd pointed out that this model supports vor- 
tices as shown in figure U2. These change color but pre- 
serve their handedness as they diffuse, and two defects 
of opposite handedness can annihilate when they meet. 
There are also chargeless excitations, where sites to ei- 
ther side of the ends of the frustrated bond have the 
same color. These can be annealed away without inter- 
acting with other defects, so they disappear exponentially 
quickly in a quench. 

Kolafa defines the charge within a region as a sum 
around a counterclockwise perimeter. 
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where m,(k) — 0, -1-1 or —I, and m{k) = k (mod 3). 
Equivalently, we can define a complex local magneti- 
zation ^i = (— l)Pe^'^**'/'^ where the three colors corre- 
spond to vectors 27r/3 apart, and (—1)^ gives a sign of 
— 1 and -|-1 on odd and even lattice sites respectively 
|ll| . This makes regions with antiferromagnetic order, 
say with color 1 on one sublattice and colors 2 and 3 
on the other, relatively uniform. Then Q is the winding 
number of ^^ around the origin. 
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FIG. 12. Counterclockwise and clockwise delects maintain 
their handedness as they diffuse. Pairs ol opposite handed- 
ness can annihilate or turn into chargeless excitations, which 
diffuse until one end is surrounded by only one other color, 
at which point it can disappear by changing to the third one. 

In figures n3 and 14, we show the time average of ^ for 
a configuration with two defects. As with the triangu- 
lar antiferromagnet, we can clearly see how each vortex 
produces a field around it which falls off with distance. 

Of all these models, this one has the easiest argument 
for a Coulomb force. As before, we can define a field E 
perpendicular to each bond which will tend to point to- 
wards clockwise defects and away from counterclockwise 
ones. We do this by orienting it so that the higher color in 
the cyclic ordering 3 > 2 > 1 > 3 is on its left. As figure 
Wb shows, after we flip one site of a positive (counter- 
clockwise) defect, the displacement of its midpoint from 
that site is just one-half the field on the edge extending 
the frustrated bond. Since both sites are equally likely, 
the average movement of the defect's midpoint is propor- 
tional to E averaged over the two bonds extending from 
the defect, 
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Since we can rewrite Kolafa's formula as 
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we again have a first-order \/r force if E is sufficiently un- 
correlated around a large perimeter. Assuming isotropy 
gives 
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FIG. 13. The local magnetization around a pair of vortices 
in the g = 3 Potts antiferromagnet. This picture was obtained 
by relaxing a 32 x 32 lattice until only two defects remained, 
and then holding them fixed for 10^ updates per site while we 
averaged the magnetization at each site. 



in the local mean-field approximation used in the previ- 
ous two models where (1/2) 'Y^E = E. 

Lenard [p3 has pointed out that this field maps the 
q — i square Potts antiferromagnet onto the six-vertex 
ice model. When defects are present, this mapping breaks 
down as shown in figure |l^. The frustrated bond has 
three outgoing or incoming arrows at each end. Thus 
defects act like monopolcs as in the previous section, but 
this time composed of a bound pair of vertices with an 
unoriented bond between them. Nijs et al. ^^ defined a 
27-vertex model that includes bonds and vertices of this 
type and studied its scaling properties. 

As for the triangular antiferromagnet, the Potts anti- 
ferromagnet can also be mapped to a solid-on-solid model 
where each site is given a height hi such that hi = Si 
(mod 3) and hi differs by ±1 between neighbors. Defects 
then correspond to a bound pair of screw dislocations 
with a total Burgers vector of ±6. 

Finally, we note that Bakaev and Kabanovich J2| have 
discussed the motion of a different kind of defect in the 
q — 2t Potts antiferromagnet, a hole with an undefined 
color. 
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FIG. 14. A scatterplot of the time-averaged local magneti- 
zation in figure [13. 
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FIG. 15. If we define a field £ on each bond so that the 
higher color in the cyclic ordering is on its left, the movement 
of a positive defect is proportional to E averaged over the two 
dashed edges. 
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FIG. 16. Mapping the q ~ S Potts antiferromagnet onto 
the six-vertex ice model and the loop-covering model. In the 
former, the defect corresponds to a monopole, an undefined 
edge with six arrows emerging from it. 




FIG. 17. A plot of tn{t) vs. log^Q t for the triangular anti- 
ferromagnet. The straight line suggests that n(t) ex (lni)/i 
for large t, and extends about a decade and a half. The data 
was taken by averaging 100 trials of 10''' updates per site on 
a 4096 X 4096 lattice, at which time ~ 75 defects remained. 



V. RELAXATION AND DIRECT 
MEASUREMENTS OF THE FORCE 

Yurke et al. [p^ discuss the relaxation dynamics of 
the XY model, in which vortices of opposite type are 
attracted with a Coulomb force. Assuming a viscous dy- 
namics V = TF where the mobility T is inversely propor- 
tional to Inr, they use simple scaling arguments to show 
that the defect density p obeys 



1 dp _ C 
p2 dt \ii{p/pc) 
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where pc is a core density and C is a constant. Note 
that the left-hand side of equation || is constant for the 
mean- field behavior p ex t^^. The leading behavior of p 
as a function of time is then 
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modifying the asymptotic behavior p ^ t~^ with a log- 
arithmic correction. They confirm the existence of this 
correction through numerical experiments. 



We performed similar experiments for these three mod- 
els. In each case, we quenched the system from T = oo [a, 
random initial state) to T = 0, and measured the num- 



ber of defects n{t) as a function of time. In figure 17 we 
graph tn{t) against logi for the triangular antiferromag- 
net, averaged over 100 runs of 10^ updates per site each 
on a 4096 x 4096 lattice. While logarithmic corrections 
are difficult to establish numerically, n{t) seems to behave 
as (lni)/i over one and a half decades in t. In addition, 
in figure |l^ we use the same data to graph p^{dp/dt)~^ 
and fit it to C~^ \n.{p/pc) as per equation |[ We obtain a 
good fit over a decade and half in p, with a core density 
Pc = 0.02. Results from the other two models are similar. 

However, this turns out not to be a good test for a 
Coulomb force. We have performed a continuous-time 
free diffusion experiment, in which we use a random state 
of the triangular antiferromagnet to set up the initial 
distribution of defects. We then diffuse the defects with 
random walks, moving them up, down, left and right 
with equal probabilities, and annihilate pairs of opposite 
charge when they meet. 

Since the total charge in an region of size I is a sum of 
0{l) surface terms, rather than 0{P) random variables, 
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FIG. 18. A plot of p^ / {dp / dt) vs. logj^g p for the triangular 
antiferromagnet, using the same data as in figure WA. The 
derivative at each time t was defined by a linear fit to 101 data 
points centered around t. The fit is to the form in equation 
|9[ If the mean-field behavior p oc t~^ held, the graph would 
be a horizontal line. 

the total charge fluctuates as 0{l^/'^) rather than 0{l). 
Thus defects of opposite charge are well mixed with each 
other in the initial condition, giving a decay close to the 
mean-field behavior p(t) oc t~^. With uncorrelated ini- 
tial conditions, defects of like type clump into domains, 
giving a t~^/^ decay p3|,p|. 

In figure |l^ we graph tn{t) for a 4096 x 4096 lattice over 
a continuous time of 8 • 10^ updates per site. The data 
appears to have the same asymptotic behavior of (lnt)/i, 
and fits of p^ /{dp/dt) to Inp are at least as convincing 
as for our three models. In other words, the asymptotic 
behavior of p is more a consequence of the correlations 
in the defects' initial positions than of the forces between 
them. This is presumably because a l/r force causes the 
length scale of the system to grow asymptotically as t^/'^, 
no faster than diffusion would anyway. At short times, 
however, p decreases faster in our three models than in 
the free diffusion experiment, indicating that attractive 
forces play a role early on when defects are relatively 
close together. 

Even if the behavior of p at long times doesn't confirm 
the existence of a Coulomb force, we can argue that it 
does show that the force between defects does not fall 
off more slowly than l/r. Generalizing the argument of 
IpTj , if the force between two defects goes as r~" for some 
a < \ and the mobility is roughly constant, the typical 
velocity is 

where ^ is the length scale of the system. Then ^ oc 
^1/(1+0)^ and the defect density is 




FIG. 19. A plot of tn{t) vs. logj^p t in a continuous-time free 
diffusion experiment, where the triangular antiferromagnet 
was used to set up the initial distribution of defects. The 
straight line suggests that n{t) oc (lnf)/f, and extends three 
decades. The data was taken by averaging 100 trials of 8 ■ 10^ 
updates per site each on a 4096 x 4096 lattice until ~ 100 
defects remained. 

p oc r' oc i-2/(l+") 

If ct < 1, then, p would decay faster than t^^ . Since we do 
not observe this, we claim that the field lines spread out, 
and are not concentrated between defects. (In another 
paper, in progress, we argue that field lines do bunch 
into tubes when a ferromagnetic next-nearest-neighbor 
interaction is added.) 

In the case of the g = 3 square Potts antiferromagnet, 
we have also measured the force between defects directly. 
In an effort to sample the set of configurations with two 
defects a particular distance away from each other, we set 
up an initial condition with two vertical defects of oppo- 
site type as in figure ^, and then let the lattice evolve 
while keeping these defects fixed. We did this on lattice 
sizes of 32, 64, 128 and 256, averaging over 10^ updates 
per site in each case after an initial equilibration of 10'^ 
updates per site. Letting the interdefect distance range 
up to 1/3 the lattice size, we perform a least-squares fit 
to the form 



(Af) 



A 



ra 



where rg is a core radius. This fit seems to work fairly 
well. 

Our values for A and tq for various lattice sizes L are 



L A 



To 



32 1.142 2.573 

64 0.978 1.956 

128 0.906 1.652 

256 0.869 1.494 



These show considerable finite-size effects. Since the 
model is critical at T = 0, we assume that A converges to 




(log-log), L=256 



FIG. 20. An initial condition with two vertical defects an 
even distance away from each other. The six checkerboard 
phases cycle around the defects, meeting at angles of 7r/3. 



its exact value as a power law in the lattice size. A least- 
squares fit to the form A{L) = A+CL" gives an exponent 
of a = —1.137 and an extrapolated value of A ^ 0.843. 
This is larger than the local mean-field approximation 
3/27r ~ 0.477 given in the previous section. We claim 
below that the exact value is A = 3/4, in which case 
our agreement is reasonable but not particularly good. 
However, this extrapolation is based on only four differ- 
ent lattice sizes, so it should be easy to improve it with 
further numerical work. 

We suspect that A is larger on smaller lattices because 
field lines have less room to spread out in. We would 
like to analyze this in terms of image charges, but it is 
somewhat unclear how to do this. In fact, even when the 
interdefect distance is half the lattice size, our prediction 
of A/{r + ro) for the force is only 25% off, even though 
by symmetry the image charges should cancel and make 
the force zero. The reason for this is that single spin- 
flip dynamics is not quite ergodic at zero temperature 
when the defects are held fixed; the boundaries between 
checkerboard domains in figure ^ can move and join with 
droplet boundaries, but not cross. Therefore, the topol- 
ogy of the field lines stays the same, with five connecting 
the defects directly and only one going around the lattice 
the other way. In a sense, this may be a happy accident, 
allowing us to observe the force at larger distances than 
we would be able to on finite lattices if image charges 
were fully felt. It also means that the measured value of 
A may depend on the topology of the initial condition; 
however, experiments on lattices of size 64 and 128 where 
all six field lines connect the two defects directly give val- 
ues of A and vq differing from those above by less than 
0.4%. 

We also tried to use the Wang-Swendsen-Kotecky clus- 



FIG. 21. The force between two defects as a function of 
their distance r , This data was obtained on a lattice of size 256 
by starting with the initial condition shown in figure EO and 
averaging the force over lO'^ updates per site after an initial 
equilibration of 10^ updates per site. During these updates, 
the defects are kept fixed, and (Ax) is defined according to 
which way they would move if we allowed them to. The form 
A/{r + ro) fits the data within the expected error over most 
of the range. 



ter algorithm |34| to sample the set of configurations 
with a pair of defects at particular places. Unfortu- 
nately, while cluster-fiipping moves leave defects in the 
same place, they sometimes replace two defects of op- 
posite charge with two chargeless ones. Finding an al- 
gorithm more efficient than single spin-flip dynamics to 
sample this set is still an open question. 



VI. ENTROPIC COULOMB FORCES 

In the XY model and other two-dimensional models 
with continuous degrees of freedom, the energy t/ of a 
free defect is typically proportional to In r where r is the 
interdefect distance [£1| . Then the energy gradient gives 
an attractive Coulomb force between defects of opposite 
type, 

dU 1 

F = -— - ex — 

or r 

In these discrete models, on the other hand, the energy 
of a pair of defects is 2 regardless of how far apart they 
are. Therefore, the force must be driven by a gradient 
in the entropy, or equivalently, by a gradient in the free 
energy at a nonzero effective temperature. 

To show how an entropy gradient can drive a first- 
order force, suppose that we group the set of spin config- 
urations with two defects into a set of macrostates S(r), 
one for each interdefect distance r. As the defects move 
towards and away from each other, we can describe the 
system as a biased random walk 



• • • ^ E(r - 1) ^ E(r) ^ S(r + 1) ^ • • • 

If we assume that this walk is a Markov process which 
has not yet had time to hit the absorbing state S(0) 
where the defects annihilate, and if during this tran- 
sient all microstates with two defects are equally likely, 
then the ratio of transition probabilities between neigh- 
boring macrostates must be the ratio of the number of 
microstates in them, 



P{r + 1 



n{r) 



where fi{r) is the number of microstates with defects r 
apart. The average motion is then 



(Ar) = P{r 



2r 



n{r 



r + 1)- P{r 



1 



n{r + l) + nir) 



r) 

1 dn{r) 
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where 



(1/2) (P(r ^ r + 1) + P{r + 1 ^ r)) 



can be thought of as a mobility, and 



F = 



dSjr) 
dr 



(10) 



(11) 



is the entropic force where S" = lnf2. 

To get a force proportional to l/r, we need to show 
that the presence of a defect decreases the entropy by an 
amount proportional to Inr. For the dimer model on the 
square lattice, Fisher and Stephenson ||l^ used Pfafhans 
and Toeplitz determinants to show that the number of 
configurations with two 'holes' a distance r apart is re- 
duced by v4r~^/^, so AS* = —(1/2) Inr. loffe and Larkin 
Jl7[ pointed out that this leads to a Coulomb force. In 
addition, Zeng et al. p9] performed numerical calcula- 
tions of the cost of dislocation pairs in a randomly pinned 
fully-packed loop model, and found a cost proportional 
to Inr. We will give several arguments for a similar term 
in the models discussed here. 

One simple counting argument goes as follows. Sup- 
pose we are trying to extend a spin configuration out- 
ward from a square region of the lattice by adding an 
additional layer of I sites around its perimeter. For the 
Potts model, for instance, successive sites have colors Si 
differing by ±1, and the sum in equation is the total 
displacement of a walk of length I. The number of such 
walks is 

I 
V2-6( 

where Q is the charge inside the perimeter. The larger 
Q is, the greater the constraint on the walks, and the 
lower the entropy will be. If we ignore the interaction 
between successive layers, which of course we can't, then 



the entropy of the set of states surrounding a charge Q 
(we take Boltzmann's constant to be 1) is 



^Q^lnH 



/=o 



El- 



I 

1/2 ^QQ 

I 
1/2 -QQ 



= 5o-E(ln(j2)-l^^(//2-6Q 



S,~{QQfJ2\ 



1=0 

for some constant A. We cut off our sum at some max- 
imum perimeter proportional to the interdefect distance 
r, outside which the charge Q is cancelled by defects of 
the opposite type. 

While the assumption of independent layers is hugely 
wrong (for instance, it gives a ground state entropy per 
site of So = Sq/N = 2) we find that the presence of a 
charge Q reduces the entropy by an amount AS' propor- 
tional to —Q^ In r. We can think of this as a contribution 
to an effective free energy. 



G 



-AS ^ AQ^lnr 



This matches nicely with the energy of a dislocation with 
Burgers vector B, 



U ^CB^lnr 



(12) 



where C is a constant dependent on the lattice spacing, 
the rigidity modulus, and Poisson's ratio [ p4| . 

As another approach, Kotecky |22] points out that the 
entropy of the Potts model can be related to the energy 
of a ferromagnetic Ising model at a nonzero effective tem- 
perature, since the entropy of one sublattice is higher if 
sites on the other sublattice are the same. An odd site 
(say) has two choices of color if its four neighbors have 
the same color, or equivalcntly if the field E points clock- 
wise or counterclockwise around all four bonds, and only 
one choice otherwise. Recall that the field E determines 
whether the color changes by -fl or —1 between neigh- 
bors. If we define the probability of the color increasing 
in the x and y directions as Px and py respectively, we 
have 

P, = (1 - Ey)/2 

Py^{l + Ex)/2 

where E^ and Ey are the components of the average field. 
If E is slowly varying, and if the colors of the four neigh- 
bors are independent, which they aren't, the probability 
of all four having the same color is 



10 



2p,{l-p.)Pyil-Py) = -{1-El)il-Eyf « -il~E') 



8 



8 



for small fields E <^ 1. The average entropy per site is 



In 2 



In 2 



-il-E^) = so~-E^ 

where sq is a (badly underestimated) ground state en- 
tropy of (1/8) In 2 per site. The effective free energy is 
then increased by 



G=i^ / E^dxdy 
8 



(13) 



giving an energy density proportional to the square of 
the field, just as for electromagnctism. 

To make this more precise, we can use the height rep- 
resentation for these models |28|, in which there is an 
effective free energy 



G^-K I \Vh\'^dxdy 



Note that V/i is simply E rotated counterclockwise by 
7r/2, so |V/ip = E^ and this has the same form as equa- 
tion |l^. A screw dislocation with Burgers vector B has 
a field around it 



IV/^I ^\E\^^^ 



B 

2T:r 



so the free energy of a defect integrated from a short- 
distance cutoff ro (roughly the lattice spacing) to an in- 
terdefect distance r is 



G 



KB^ 

47r 



(In?' — Inro) 



which again matches the form in equation y_2[ The force 
between two screw dislocations with Burgers vector B 
and B' a distance r apart is then 



F ■ 



dG 
dr 



K- 



BB' 



Combining this with equation 11 gives 

(Ar) = 



(14) 



(15) 



where B is the magnitude of the Burgers vector for a 
defect of charge ±1. 

For the q — 2> Potts antiferromagnet, we can calculate 
the coefficient A = TKB^/n of equation na exactly. Park 
and Widom BO] , using the exact solution of the six- vertex 
ice model, showed that the free energy of an interface of 
width L across a height difference of Ah = ±2 is 27r/6L, 
and Burton and Henley [nO[ point out that setting this 
equal to {l/2)K{Ah/L)'^L gives K = n/6. 

To calculate F, we need to take into account the fact 
that a diffusing defect alternates between horizontal and 



vertical bonds. If we assume for simplicity that E is in the 
x-direction, a vertical defect which flips to a neighboring 
horizontal bond has an equal probability of flipping back 
to its original position or moving to the vertical bond 
one site away. If P^ and P^ are the probabilities that 
it will move to the horizontal bond to its right and left 
respectively, after two time-steps we have P{r — > r + l) = 
P^/2 and P{r + I ^ r) = P^/2. Since P^ + P^ ^ 1, 
plugging these into equation ^ and dividing by 2 since 
F is inversely proportional to time gives F — 1/8. 

Combining these with B = 6 gives A = 3/4. Our 
experimental value of 0.841 differs from this by 12%, but 
this extrapolation was based on only four lattice sizes. 
We hope to measure A more accurately in the future. 



VII. CONCLUSION 

We have examined topological defects or vortices in 
the triangular Ising antiferromagnet, a model with two- 
and four-point interactions on the square lattice equiva- 
lent to the six-vertex ice model under a dynamics where 
elementary moves reverse the four arrows around a pla- 
quet, and the three-state Potts antiferromagnet on the 
square lattice. In each case, positive and negative de- 
fects appear to be attracted by a Coulomb force, and we 
have confirmed this through numerical experiments. For 
the Potts model, we have obtained reasonable agreement 
between the force coefficient and its predicted value. 

Both the triangular Ising antiferromagnet p3,p6[ and 
the g = 3 Potts antiferromagnet 129] are known to 
have Kosterlitz-Thouless-like phase transitions |2^ in 
the same universality class as the six-state clock model 
p9| , p6[ when a ferromagnetic next-nearest-neighbor in- 
teraction is added. In another paper, we will look at 
how this interaction, and the screening effect of particle- 
antiparticle pairs at nonzero temperatures, affects the 
forces between defects. We will also look at vortex loops 
in various three-dimensional models. 
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